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APPLICATION  OF  FLUX-CORRECTED  TRANSPORT  TO  TTCP 
JOINT  LAUNCH  BLAST  COMPUTATIONAL  EFFORT 

INTRODUCTION 

/ 

-  The  method  of  Flux-Corrected  Transport  (FCT)  vas  developed^gbout  ten 
years  ago  (Boris  and  Book,  1973;  Book  et  al. ,  1975;  Boris  and  Book,  1976l^as 
a  means  for  solving  systems  of  hyperbolic  equations  in  such  a  way  that 
physically  positive  quantities  like  mass  and  energy  density  remain  positive. 
Its  principal  application  has  been  to  compressible  fluids,  i.e.,  fluids  in 
which  some  of  the  flow  speeds  are  comparable  with  or  greater  than  the  local 
speed  of  propagation  of  waves  in  the  medium.1  Examples  include  plasmas, 
combusting  systems,  and  gas-dynamic  flows.  Applications  of  FCT  to  muzzle 

blast  problems  were  first  made  by  Townend  and  Edwards  (1982). 

^ ' 

During  the  past  decade  the  method  has  been  extended  and  generalized  in  a 
number  of  important  respects."  These  include  the  construction  of  general- 
purpose  algorithms  for  solving  fluid  equations  on  moving  nonuniform  grids  in 
a  curvilinear  coordinate  system  (Boris  1976),  strictly  multidimensional  FCT 
algorithms  (as  opposed  to  timestep-splitting  of  one-dimensional  algorithms, 
Zalesak  1979 )»  and  variants  on  these.  In  addition,  numerous  techniques  have 
been  developed  by  other  workers  which  incorporate  the  positivity -preserving 
property  in  other  ways  (Harten,  1978;  Van  Leer,  1979;  Collela  and  Woodward, 
198L).  The  weight  of  opinion  among  computational  scientists  who  deal  with 
problems  of  compressible  flow  now  strongly  affirms  the  superiority  of 
positivity -preserving  methods  (e.g.  ,  Total  Variation  Diminshing  or  TVD)  over 
conventional  ones.  The  advantages  are  improved  robustness,  flexibility, 
resolution,  and  freedom  from  nonphysical  numerical  artifacts. 

Manuscript  approved  October  9,  1985. 
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FAST2D  (Boris  1976)  is  an  FCT  code  for  solving  the  two-dimensional  fluid 
equations  in  Cartesian  or  r-z  geometry.  It  employs  timestep-splitting  to 
permit  the  use  of  JPBFCT,  a  general-purpose  one-dimensional  convective 
equation  solver  for  both  coordinate  sweeps  in  all  of  the  fluid  equations. 

To  explain  what  this  means,  let  us  describe  a  class  of  FCT  algorithms 
which  includes  JPBFCT  as  a  particular  case.  To  advance  the  one-dimensional 
continuity  equation 


If-  -  -  k  (pv) 


(i) 


one  timestep  on  a  uniform  mesh  with  a  given  flow  velocity  v,  we  suppose  the 


old  cell-centered  values  p  and  v  are  known.  A  three-point  conservative 

J  J 

o 

finite-difference  operator  acting  on  p  can  be  written  in  the  form 

J 


T  o 
P.  =  P,  -  £ 


i  ~  fcj+l/2  PJ+l/2  +  GJ -1/2  pJ-l/2  » 


(2) 


where  p  .  ^  =  or  (p  +  p  .  ,,  ).  To  obtain  a  first-order  approximation  to  Eq. 

J+l/^  2  j  J+l 

(l),  we  must  define  z the  Courant  number  Vj+1^,26t/5x,  where 

1  T 

vJ+l/2  =  2  ^vj  +  VJ+1^’  In  a  numerical  diffusion  is  then  applied  to  pj 

in  flux  form  according  to 


TD  T  /  o  o  ^  /  o  o  \ 

°J  ~  PJ  j+l/2(pj+l  "  pjJ  VJ-l/2  pj  "  pj-lJ  *  (3) 

Ilext  we  act  to  eliminate  this  excess  diffusion.  If  we  simply  applied 


another  diffusion  operation,  this  time  with  the  opposite  sign,  it  would 

TD 

cancel  the  diffusion  terms  already  present  in  p  ,  Instead  we  define  new 

J 

,  . .  .  n 

densities  p 


n  TD  c  c 

PJ  "  +1/2  *J-l/2’ 


(M 


V  V  * 
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where  the  "corrected"  fluxes  *J+l/2  differ  from  the  "raw"  fluxes  t^+1/2 


wJ+l/2  ^pj+l  -  in  tWO  ways:  we  use  fpjl  instead  of  {p^}  in  the  definition 
of  the  raw  fluxes  (this  allows  us  to  optimize  some  property  in  the  difference 
scheme,  as  will  be  seen  shortly);  and  we  correct  the  fluxes,  so  that  the 
antidiffusion  process  (which  evidently  tends  to  make  all  gradients  steeper) 


TDi 

can  enhance  no  extrema  already  present  in  {p^  },  nor  introduce  any  new  ones. 


The  simplest  formula  for  achieving  this  in  all  possible  situations  is  that 
used  in  "strong  flux  limiting": 


♦S+l/2  =  S*{maxI°’  111111  (IVl/2!’  Vl/2’  V3/2)1}  ’  (5) 

TD  TD 

where  S  -  sign  *J+1/2  and  AJ+1/2  =  Pj+1  -  Pj  . 

At  most  points  in  a  gently  varying  profile  no  correction  is  required  and 

^J+l/2  =  +1/2*  In  that  case  we  can  perform  a  Von  Neumann  analysis, 

calculating  the  complex  propagator  or  amplification  factor  A  =  A  +  iA.  = 

r  i 

Pj/p°  for  a  sinusoidal  density  profile  p°  =  exp(ijk  5x),  where  k  is  the  wave 
number.  Writing  B  =  kdx  and  e  =  v6t/<5x,  we  expand  the  amplification  in 
powers  of  8 


A  =  1  +  Ag  82  +  A^  B4  +  ... 
and  the  relative  phase  error 

R  -  tan-1  (A^/A^)  -  1  =  R2  B2  +  R4  B^  +  .... 


(6) 


(7) 


It  is  easy  to  show  that  the  condition  that  A2  vanish  is 


The  condition  that  R2  vanish  is 


1  1  2 
w  =  I  "  S  e 


(9) 


While  other  choices  of  v  and  y  are  possible,  Eqs.  (8)  and  (9)  lead  to  the 
most  useful  general-purpose  algorithm  within  this  class  and  have  been  used  in 
a  wide  variety  of  supersonic  flow  problems.  Besides  producing  small  phase 
and  amplitude  errors,  as  shown  above,  this  choice  of  the  diffusion  and 
antidiffusion  coefficients  v  and  y  treats  shocks  very  accurately,  unlike 
schemes  where  v  and  y  are  small. 

Generalizing  this  scheme  to  handle  the  momentum  and  energy  equations 
presents  no  problem  as  the  treatment  of  the  driving  terms  -Vp,  etc.,  and 
-7«pv  is  not  very  sensitive.  Extension  to  nonuniform  meshes  and 
curvilinear  geometry  is  also  straightforward.  When  fully  vectorized  and 
optimized  for  a  particular  machine  (in  the  present  case,  the  Cray-1)  it  runs 
at  about  2  ys  per  zone  per  timestep  per  equation  per  coordinate  direction. 


RESULTS 

Case  I:  Flow  From  the  Open  End  of  a  Shock  Tube 

Figure  1  reproduces  the  BRL  drawings  defining  the  geometry  of  the  shock 
tube  and  prescribing  the  computational  domain  to  be  used  (Schmidt,  1985). 
Although  it  would  be  of  interest  to  model  the  entire  length  of  the  shock  tube 
(including  the  0.673  m  driver  section  and  the  whole  3.0^1  m  driven  section), 
this  would  necessitate  a  calculation  with  a  computational  region 
approximately  ten  times  as  long  as  that  specified  for  the  case.  In  our 
calculation  we  used  a  grid  with  200  zones  in  the  radial  and  LOC  zones  in  the 
axial  direction,  with  mesh  sizes  Ar  =  Az  =  0.228  cm.  Although  temperature 
and  density  do  not  deviate  too  widely  from  standard  conditions  and  y  =  l.U 


holds  approximately,  we  used  a  real-air  equation  of  state.  Reflecting 
boundary  conditions  were  employed  at  r  =  0,  and  outflow  conditions  (valid  for 
supersonic  flow)  at  maximum  radial  and  axial  displacement.  The  conditions  at 
2=0  were  maintained  at  their  initial  values,  which  is  valid  until 
disturbances  propagating  out  of  the  driver  region  reach  the  computational 
domain. 

Timesteps  were  limited  to  0.5  times  the  value  imposed  by  the  Courant 
condition  and  averaged  about  a  microsecond.  The  calculation  ran  l6l8  cycles, 
with  dumps  being  obtained  at  times  t  =  0.025,  0.05,  0.1,  0.2,  0.5,  1.0  and 
1.5  ms.  Using  these  dumps  we  generated  velocity  vector  plots  and  contours  of 
constant  density,  pressure,  and  Mach  number.  In  addition,  time  histories  of 
pressure,  density,  and  Mach  number  were  generated  at  stations  situated  at  a 
distance  of  1.5  times  the  tube  diameter  D  =  0.152  m  from  the  center  of  the 
tube  exit  along  the  0,  30,  60,  90,  120,  130,  and  138.2  degree  rays.  These 
locations  are  indicated  by  the  dots  in  Fig.  1c.  Since  the  station  along  the 
0°  ray  was  of  particular  interest  in  connection  with  determining  the 
location  of  the  backward-facing  shock,  ten-  additional  stations  were  located 
in  the  neighborhood  of  the  station  along  this  ray. 

Figures  2-8  show  the  overall  flow  at  the  plotting  times  specified  above. 
In  each  figure  we  have  four  plots:  contours  of  (a)  pressure,  (b)  density, 
and  (c)  Mach  number;  and  (d)  velocity  vector  plots.  Note  that  by  0.5  ms 
(Fig.  6)  several  well-defined  gasdynamic  discontinuities  have  developed. 

These  include  a  disk-shaped  backward-facing  shock  which  terminates  at  its 
outer  radius  in  a  compression  fan.  The  latter  is  rolled  up  into  a  vortex 
ring.  At  a  radius  of  approximately  0.05  m  a  cylindrical  slip  surface  (which 
shows  up  as  a  discontinuity  in  density  but  not  pressure)  intersects  this 
shock,  producing  a  kink.  By  1.5  ms  (Fig.  8),  when  the  primary  shock  has 


^  w  v  V 


i  —  *•  “V  *M 


reached  the  end  of  the  mesh,  the  vortex  ring  has  expanded  out  to  about  0.1  m 
in  radius. 

Figures  9  -  16  are  station  plots,  each  showing  (a)  overpressure,  (b) 
density,  and  (c)  Mach  number  at  the  prescribed  locations.  The  flow  behind 
the  shock  is  of  course  subsonic  everywhere.  All  stations  show  a  pressure 
(and  density)  peak  at  shock  arrival  time,  followed  by  a  rarefaction  wave. 

The  stations  located  off  the  center  line  (all  except  no.  l)  reach  pressures 
less  than  ambient.  Comparison  with  digitized  traces  of  pressure  measurements 
(Figs.  17  -  33)  shows  good  agreement  between  calculated  and  recorded  values 
except  at  station  no.  1.  There  the  calculated  values  (labelled  with  "A")  are 
close  to  the  measured  ones  (labelled  with  "B")  until  about  1  ms.  Thereafter 
the  latter  peak  up  above  the  calculation  by  about  20  kPa  and  then  drop  below 
by  about  60  kPa.  We  attribute  this  discrepancy  to  slight  errors  in  defining 
the  location  of  the  backward-facing  shock,  probably  caused  by  an  increase  in 
the  effective  flow  velocity  from  the  opening  of  the  shock  tube  as  a  result  of 
the  formation  of  a  boundary  layer  there.  Rather  than  try  to  model  this 
boundary  layer,  we  introduced  ten  additional  stations  on  the  axis  between 
C.662  m  and  0.683  m  from  the  opening.  The  waveform  which  agrees  best  with 
the  observed  trace  is  that  calculated  at  station  no.  10  (Fig.  26,  y  =  0.667 
m).  Nevertheless,  the  unsteady  behavior  of  the  shock  in  the  experimental 
setup  precludes  perfect  agreement  for  these  features. 

Case  II:  10$  mm  Muzzle  Flow 

In  this  problem  the  geometry  is  complicated  by  the  addition  of  a  baffle 
located  downstream  from  the  muzzle,  and  by  the  presence  of  the  projectile  at 
early  times.  No  attempt  was  made  to  model  the  latter  (although  it  is  r.ot 


difficult  to  do)  because  it  has  no  influence  on  the  flow  near  the  nuzzle 


except  at  very  early  times.  Otherwise  the  same  boundary  conditions  were  used 
as  in  Case  I. 


«\  A  A 
.■•v  •.* ' 


The  zone  size  was  Ar  =  Az  =  0.3  cm.  For  this  problem  two  fluid  mass 
densities  were  propagated  according  to  the  continuity  equation,  using  a 
single  flow  velocity  field.  One  fluid,  air,  initially  filled  the  entire  grid 
outside  the  muzzle  and  the  solid  walls.  The  other  fluid,  representing 
detonation  products  from  the  propellant,  expands  out  from  the  muzzle  with  the 
prescribed  exit  velocity  of  1050  m/s.  Partial  pressures  were  obtained  using 
respectively  a  real-air  equation  of  state  and  a  constant-gamma  law  (y  =  1.25) 
for  the  two  fields. 

Figures  31*  -  38  shows  contours  of  (a)  pressure,  (b)  density,  and  (c) 

Mach  number,  and  (d)  velocity  vector  plots  at  the  prescribed  times  t  =  0.025, 
0.05,  0.1,  0.5,  and  1.0  ms,  respectively.  Note  that  by  0.1  ms  part  of  the 
shock  wave  has  reached  the  baffle  and  been  deflected  away.  The  central 
portion  passes  through  the  hole  and  begins  to  diffract  around  to  the  back  of 
the  baffle.  By  0.5  ms  the  shock  has  left  the  top  of  the  grid.  Thereafter 
the  solution  quickly  relaxes  to  a  steady  state,  maintained  by  the  propellant 
gases  emerging  from  the  howitzer  muzzle. 

Figures  39  -  **7  show  waveforms  of  (a)  pressure,  (b)  density,  and  (c) 

Mach  number  at  the  nine  prescribed  stations.  Stations  (l)  -  (3),  located  on 
the  front  of  the  baffle,  show  slightly  different  peak  values  at  shock  arrival 
time  and  saturate  at  slightly  different  values  at  late  times,  reflecting 
their  varying  distances  from  the  muzzle.  The  stations  in  front  of, 
alongside,  and  behind  the  muzzle  show  variations  which  can  readily  be 
interpreted  in  terms  of  geometry.  It  is  seen  that  although  the  flow  is 
deflected  by  the  baffle,  beyond  it  (at  distances  >  D)  the  baffle  has  little 
effect  on  most  features  in  either  contour  or  station  plots. 


'•V-V-'j 
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Since  the  flow  is  subsonic,  the  mesh  boundaries  should  permit  inflow  of 
the  surrounding  ambient  atmosphere  at  late  times,  and  this  should  come  to 
equilibrium  with  the  emerging  flow  of  propellent  gases.  Because  we  kept  the 
same  outflow  boundary  conditions  as  in  Case  I,  inflow  cannot  take  place,  and 
at  late  times  the  pressures  and  densities  near  the  outer  boundary  become 
unphysically  low.  For  this  reason  we  have  plotted  station  histories  only  out 


CONCLUSION 


We  have  carried  out  two  test  muzzle  blast  calculations  using  the  FCT 
code  FAST2D.  Numerous  simplifications  have  been  introduced  to  minimize  the 
programming  effort  required,  and  considerable  improvement  could  be  made  in 
our  treatment  of  (i)  viscous  boundary  layers  (in  Case  I);  (ii)  inflow  and 
projectile  boundary  conditions  (in  Case  II);  and  (iii)  the  propellant  gas 
equation  of  state.  Nevertheless,  the  agreement  with  the  data  in  Case  I  and 
the  physical  plausibility  and  internal  consistency  of  both  cases  strongly 
support  the  validity  of  this  modeling  approach  and  suggest  that  FCT  hydro¬ 
codes  can  make  a  valuable  contribution  to  the  muzzle  blast  modeling  effort. 

Although  the  calculations  described  here  were  carried  out  on  a  Cray-1 
supercomputer,  it  is  worth  noting  that  the  computing  requirements  for  FAST2D 
are  actually  very  modest.  The  Laboratory  for  Computational  Physics  has  as¬ 
sembled  a  Graphical  and  Array  Processing  System  (GAPS)  consisting  of  a  VAX 
11/780,  two  FPS  5305  array  processors  working  in  tandem  (up  to  10  APs  are 
possible),  an  Aptec  DPS  2h00  I/O  computer,  and  a  Tektronix  U115B  high-speed 
graphics  device.  This  system  is  running  a  comparable  two-dimensional  hydro¬ 
code,  programmed  entirely  in  FORTRAN,  at  speeds  of  ~  0.2  ms  per  mesh  point 
per  timestep,  in  either  batch  or  interactive  mode.  In  the  latter  it  produces 
color  plots  of  the  evolving  simulation  approximately  every  5  seconds  with 
DICCMED  hard  copy  capability.  This  is  about  1/8  the  speed  of  a  single 
processor  Cray  X-MP.  A  sample  frame  from  a  muzzle  blast  calculation  done 
using  the  RFM  hydrocode  on  the  GAPS  is  reproduced  as  Fig.  b8.  It  is 
important  to  stress  that  such  a  system  is  small  and  can  be  maintained  in  a 
secure  self-contained  facility,  where  it  can  be  dedicated  entirely  to  muzzle 
blast  or  similar  energetic  gasdynamic  calculations  if  so  desired.  The  cost 
of  the  basic  system  is  ~  $?2CK. 
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